Ethanol

standard_theme <- theme_bw() +
  theme(plot.title = element_text(size = 30, hjust = 0.5),
    plot.subtitle = element_text(size = 25, hjust = 0.5),
    plot.caption = element_text(size = 20),
    axis.title = element_text(size = 25),
    axis.text.y = element_text(size = 20),
    axis.text.x = element_text(size = 20),
    legend.title = element_text(size = 20),
    legend.text = element_text(size = 22),
    strip.text = element_text(size = 22))
theme_set(standard_theme)

set_flextable_defaults(table.layout = "autofit")

stat <- list(
        "Mean" = ~mean(., na.rm = TRUE),
        "Std_dev"   = ~sd(., na.rm = TRUE)
      )
stat_ph <- list(
        "Среднее" = ~mean(., na.rm = TRUE) %>% round(4),
        "Медиана" = ~median(., na.rm = TRUE)%>% round(4),
        "Стандартное отклонение"   = ~sd(., na.rm = TRUE) %>% round(4),
        "Коэффициент вариации (CV), %" = ~round((sd(., na.rm = TRUE))/(mean(., na.rm = TRUE)) * 100, 2)
        )

1 Чтение данных

ethanol <- read_delim("./data/raw/Jones-1984.csv")

## Rows: 480 Columns: 9
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (5): study_id, route, food, site, sex
## dbl (4): subject_id, time_min, conc_g_per_L, dose_g_per_kg
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

ethanol <- ethanol %>% 
  mutate(across(where(is.character), as.factor))
summary(ethanol)

##        study_id     subject_id       time_min    conc_g_per_L    dose_g_per_kg 
##  Jones_1984:480   Min.   : 1.00   Min.   :  0   Min.   :0.0000   Min.   :0.68  
##                   1st Qu.:12.75   1st Qu.: 60   1st Qu.:0.2575   1st Qu.:0.68  
##                   Median :24.50   Median :150   Median :0.5300   Median :0.68  
##                   Mean   :24.50   Mean   :180   Mean   :0.4966   Mean   :0.68  
##                   3rd Qu.:36.25   3rd Qu.:300   3rd Qu.:0.7300   3rd Qu.:0.68  
##                   Max.   :48.00   Max.   :420   Max.   :1.3000   Max.   :0.68  
##  route        food            site       sex     
##  PO:480   fasted:480   capillary:480   male:480  
##                                                  
##                                                  
##                                                  
##                                                  
## 

str(ethanol)

## tibble [480 × 9] (S3: tbl_df/tbl/data.frame)
##  $ study_id     : Factor w/ 1 level "Jones_1984": 1 1 1 1 1 1 1 1 1 1 ...
##  $ subject_id   : num [1:480] 1 1 1 1 1 1 1 1 1 1 ...
##  $ time_min     : num [1:480] 0 30 60 90 120 180 240 300 360 420 ...
##  $ conc_g_per_L : num [1:480] 0 0.94 0.81 0.73 0.65 0.56 0.44 0.34 0.17 0.01 ...
##  $ dose_g_per_kg: num [1:480] 0.68 0.68 0.68 0.68 0.68 0.68 0.68 0.68 0.68 0.68 ...
##  $ route        : Factor w/ 1 level "PO": 1 1 1 1 1 1 1 1 1 1 ...
##  $ food         : Factor w/ 1 level "fasted": 1 1 1 1 1 1 1 1 1 1 ...
##  $ site         : Factor w/ 1 level "capillary": 1 1 1 1 1 1 1 1 1 1 ...
##  $ sex          : Factor w/ 1 level "male": 1 1 1 1 1 1 1 1 1 1 ...

2 Визуализация

2.1 Индивидуальные кривые концентрации-время

ggplot(ethanol, aes(x = time_min, y = conc_g_per_L)) +
  geom_line(aes(group = subject_id)) +
  geom_smooth(method = "gam") +
  scale_y_continuous(breaks = seq(0, 1.4, 0.1)) +
  scale_x_continuous(breaks = seq(0, 420, 30)) +
  # coord_cartesian(xlim = c(60, 420)) +
  labs(title = "Спагетти-плот - кривые концентрации время",
       subtitle = "Для демонстрации тренда применён метод GAM",
       caption = "Jones-1984",
       x = "Время, мин",
       y = "Концентрация, г/л")

## `geom_smooth()` using formula = 'y ~ s(x, bs = "cs")'

## Общая кривая (среднее плюс стандартное отклонение)

ethanol %>% 
  mutate(time_min = time_min %>% as.character) %>% 
  group_by(time_min) %>% 
  summarise(across(conc_g_per_L, stat, .names = "{.fn}")) %>% 
  ungroup() %>% 
  mutate(time_min = time_min %>% as.numeric) %>% 
  ggplot(aes(y = Mean, x = time_min)) +
  geom_point() +
  geom_line(aes(group = "")) +
  geom_errorbar(aes(ymin = Mean - Std_dev,
                    ymax = Mean + Std_dev),
                width = 20) +
   scale_y_continuous(breaks = seq(0, 1.4, 0.1)) +
  scale_x_continuous(breaks = seq(0, 420, 30)) +
  labs(title = "Среднее ± стандартное отклонение\nдля каждой временной точки",
       caption = "Jones-1984",
       x = "Время, мин",
       y = "Концентрация, г/л")

ggplot(ethanol, aes(x = time_min, y = conc_g_per_L)) +
  geom_smooth(method = "lm", aes(group = subject_id), se = FALSE) +
  scale_y_continuous(breaks = seq(0, 1.4, 0.1)) +
  scale_x_continuous(breaks = seq(0, 420, 30)) + 
  labs(title = "Линейная регрессия концентрации для каждого субъекта",,
       caption = "Jones-1984",
       x = "Время, мин",
       y = "Концентрация, г/л")

## `geom_smooth()` using formula = 'y ~ x'

3 Расчёт основных фармакокинетических показателей

#k - это насколько мы прибавляем к индексу tmax
pharma_fun <- function(conc, t, dose, k = 2) {
  # Моделирование для подсчёта beta, C0, Vd
  t_sample <- {{t}}[(first(which.max({{conc}})) + k):length({{t}})]
  conc_sample <- conc[(first(which.max({{conc}})) + k):length({{conc}})]
  dose <- first({{dose}})
  model <- lm({{conc_sample}} ~ {{t_sample}})
  Beta <- -as.numeric(model$coefficients[2])
  C_0 <- as.numeric(model$coefficients[1])
  V_d <- dose / C_0
  
  # Подсчёт AUC0-t
  
  # if (length({{t}}) < 2) AUC <- NA
  # if (anyNA({{t}}) | anyNA({{conc}})) AUC <- NaN
  t0 = {{t}}[-length({{t}})]
  t1 = {{t}}[-1]
  c0 = {{conc}}[-length({{conc}})]
  c1 = {{conc}}[-1]
  AUC = sum(((c0 + c1)/2)*(t1 - t0))
  
  # Подсчёт Cmax
  cmax <- max(conc)
  
  # Подсчёт Tmax
  tmax <- t[first(which.max(conc))]
  return(list(Beta = Beta,
              C_0 = C_0,
              V_d = V_d,
              AUC = AUC,
              Cmax = cmax,
              Tmax = tmax))
}

ethanol %>% 
  arrange(time_min) %>% 
  nest(.by = subject_id, .key = "subject_data") %>% 
  mutate(pharma = purrr::map(subject_data,
                            ~pharma_fun(.$conc_g_per_L,
                                        .$time_min,
                                        .$dose_g_per_kg))) %>% 
  unnest_wider(col = pharma) %>% 
  select(-subject_data) -> pharmacokinetics_data

4 Описательные статистики для фармакокинетических показателей

pharmacokinetics_data <- pharmacokinetics_data %>% 
  mutate(Beta = Beta * 60) %>%  # перевожу в г/л в час
  rename(`β, г/л в час` = Beta,
         `C0, г/л` = C_0,
         `Vd, л/кг` = V_d,
         `AUC0-t, г*мин/л` = AUC,
         `Cmax, г/л` = Cmax,
         `Tmax, мин` = Tmax)
pharmacokinetics_data %>% 
  pivot_longer(!subject_id,
               names_to = "Показатель",
               values_to = "Значение") %>% 
  group_by(`Показатель`) %>% 
  summarise(across(`Значение`, stat_ph, .names = "{.fn}")) %>% 
  ungroup() %>% 
  flextable() %>%
  theme_box %>% 
  fontsize(size=9) %>% 
  add_header_row(values  = "Описательные статистики фармакокинетических показателей", colwidths = 5)

Описательные статистики фармакокинетических показателей
Показатель	Среднее	Медиана	Стандартное отклонение	Коэффициент вариации (CV), %
AUC0-t, г*мин/л	212.3156	210.4500	25.1503	11.85
C0, г/л	0.9886	0.9806	0.0747	7.56
Cmax, г/л	0.9250	0.9300	0.1550	16.76
Tmax, мин	54.3750	60.0000	28.1263	51.73
Vd, л/кг	0.6916	0.6935	0.0512	7.41
β, г/л в час	0.1238	0.1243	0.0123	9.94

5 Построение графиков

hist_fn <- function(data, col) {
  
  ggplot(data, aes(x = {{col}})) +
  geom_boxplot(
    alpha = 0.5,
  # width = width,
    linewidth = 1.5,
    outlier.size = 2
  ) +
  geom_histogram(
    fill = "midnightblue",
    alpha = 0.5,
    colour = "black",
    bins = 25
  ) +
    geom_vline(aes(xintercept = mean({{col}}, na.rm = T)),
               color = "darkred",
               linetype = 2,
               linewidth = 2) +
  labs(# title = rlang::englue("{{col}}"),
       y = "Количество, n",
       x = rlang::englue("{{col}}")) +
  theme(plot.title = element_text(hjust = 0.5),
        axis.title = element_text(size = 25),
        axis.text.y = element_text(size = 23),
        axis.text.x = element_text(size = 23),
        legend.title = element_text(size = 23),
        legend.text = element_text(size = 23))
}

col_list <- pharmacokinetics_data %>%
  select(-subject_id) %>%
  names() %>%
  syms()

hist_plots <- purrr::map(col_list,
    ~ hist_fn(pharmacokinetics_data, !!.x))

wrap_plots(hist_plots, ncol = 2) +
  plot_annotation(title = "Гистограммы и боксплоты фарм. показателей",
                  theme = standard_theme)

6 Интерпретация результатов

1. Фармакокинетические кривые в целом лежат достаточно кучно, имеют очень похожие профили;

2. В основном \(Cmax\) достигалась на 30 минуте, всего представлены 4 различных \(Tmax\);

3. Похоже, что распределения константы элиминации и площади под кривой бимодальны;

4. Средняя \(Cmax\) составила 0.93 ⼠ 0.16 г/л;

5. Не смотря на гомогенность выборки по основным параметрам у \(Tmax\) значительный коэффициент вариации (около 50%), возможно, не учтён какой-либо фактор;

6. Мне кажется, что полученное распределения \(C_0\) и \(β\) имеют достаточно низкую вариабельность.